{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "18215fa2",
   "metadata": {},
   "source": [
    "# Tutorial 03: weak imposition of Dirichlet BCs by a Lagrange multiplier (interface problem)\n",
    "\n",
    "In this tutorial we solve the problem\n",
    "\n",
    "$$\\begin{cases}\n",
    "-\\Delta u = f, & \\text{in } \\Omega,\\\\\n",
    " u   = g, & \\text{on } \\partial\\Omega,\n",
    "\\end{cases}$$\n",
    "\n",
    "where $f=1$, $g=0$ and $\\Omega$ is a ball in 2D, using a domain decomposition approach for $\\Omega = \\Omega_1 \\cup \\Omega_2$, and introducing a lagrange multiplier to handle the continuity of the solution across\n",
    "the interface $\\Gamma$ between $\\Omega_1$ and $\\Omega_2$.\n",
    "\n",
    "The resulting weak formulation is:\n",
    "$$\n",
    "\\text{find } u_1, u_2, \\lambda \\in V(\\Omega_1) \\times V(\\Omega_2) \\times E(\\Gamma) \n",
    "$$\n",
    "s.t.\n",
    "$$\n",
    "\\int_{\\Omega_1} \\nabla u_1 \\cdot \\nabla v_1 dx +\n",
    "\\int_{\\Omega_2} \\nabla u_2 \\cdot \\nabla v_2 dx +\n",
    "\\int_{\\Gamma} \\lambda (v_1 - v_2) ds +\n",
    "\\int_{\\Gamma} \\eta  (u_1 - u_2) ds = \n",
    "\\int_{\\Omega_1} f v_1 dx +\n",
    "\\int_{\\Omega_2} f v_2 dx,\n",
    "\\qquad \\forall (v_1, v_2, \\eta) \\in V(\\Omega_1) \\times V(\\Omega_2) \\times E(\\Gamma).\n",
    "$$\n",
    "\n",
    "Equivalenty this equation can be written as a system of equations\n",
    "$$\n",
    "\\begin{align*}\n",
    "&\\int_{\\Omega_1} \\nabla u_1 \\cdot \\nabla v_1 dx & & & &+ \\int_\\Gamma \\lambda v_1 ds & &= \\int_{\\Omega_1} f v_1 ds & &\\forall v_1 \\in V(\\Omega_1) \\\\\n",
    "& & &\\int_{\\Omega_2} \\nabla u_2 \\cdot \\nabla v_2 dx & &- \\int_\\Gamma \\lambda v_2 ds & &= \\int_{\\Omega_2} f v_2 ds & &\\forall v_2 \\in V(\\Omega_2) \\\\\n",
    "&\\int_\\Gamma \\eta u_1 ds & &- \\int_\\Gamma \\eta  u_2 ds & & & &= 0 & &\\forall \\eta \\in E(\\Gamma).\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Boundary conditions on $\\partial\\Omega$ are embedded in $V(\\Omega_i) \\subset H^1(\\Omega_i)$, $i = 1, 2$, and $E(\\Gamma) \\subset L^2(\\Gamma)$.\n",
    "\n",
    "This example is a prototypical case of problems containing interface restricted variables (the Lagrange multiplier, in this case)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfd8e555",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfinx.fem\n",
    "import dolfinx.fem.petsc\n",
    "import dolfinx.io\n",
    "import gmsh\n",
    "import mpi4py.MPI\n",
    "import numpy as np\n",
    "import petsc4py.PETSc\n",
    "import ufl\n",
    "import viskex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd776df7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiphenicsx.fem\n",
    "import multiphenicsx.fem.petsc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41282b79",
   "metadata": {},
   "source": [
    "### Geometrical parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df239c5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 3\n",
    "mesh_size = 1. / 4."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aa03a00",
   "metadata": {},
   "source": [
    "### Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e34f01c",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.initialize()\n",
    "gmsh.model.add(\"mesh\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1751133",
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)\n",
    "p1 = gmsh.model.geo.addPoint(0.0, +r, 0.0, mesh_size)\n",
    "p2 = gmsh.model.geo.addPoint(0.0, -r, 0.0, mesh_size)\n",
    "c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)\n",
    "c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)\n",
    "l0 = gmsh.model.geo.addLine(p2, p1)\n",
    "line_loop_left = gmsh.model.geo.addCurveLoop([c0, l0])\n",
    "line_loop_right = gmsh.model.geo.addCurveLoop([c1, -l0])\n",
    "semicircle_left = gmsh.model.geo.addPlaneSurface([line_loop_left])\n",
    "semicircle_right = gmsh.model.geo.addPlaneSurface([line_loop_right])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "970792ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.model.geo.synchronize()\n",
    "gmsh.model.addPhysicalGroup(1, [c0, c1], 1)\n",
    "gmsh.model.addPhysicalGroup(1, [l0], 2)\n",
    "gmsh.model.addPhysicalGroup(2, [semicircle_left], 1)\n",
    "gmsh.model.addPhysicalGroup(2, [semicircle_right], 2)\n",
    "gmsh.model.mesh.generate(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "950bee0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet)\n",
    "mesh, subdomains, boundaries_and_interfaces = dolfinx.io.gmshio.model_to_mesh(\n",
    "    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2, partitioner=partitioner)\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0caef08d-fd6f-48e7-9ec6-74247edbc244",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create connectivities required by the rest of the code\n",
    "mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)\n",
    "mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)\n",
    "mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d38a7eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "cells_Omega1 = subdomains.indices[subdomains.values == 1]\n",
    "cells_Omega2 = subdomains.indices[subdomains.values == 2]\n",
    "facets_partial_Omega = boundaries_and_interfaces.indices[boundaries_and_interfaces.values == 1]\n",
    "facets_Gamma = boundaries_and_interfaces.indices[boundaries_and_interfaces.values == 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a82d6ea1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define associated measures\n",
    "dx = ufl.Measure(\"dx\", subdomain_data=subdomains)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f65f333",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_mesh(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbcbe1e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_mesh_tags(mesh, subdomains, \"subdomains\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35e6b52c",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_mesh_tags(mesh, boundaries_and_interfaces, \"boundaries and interfaces\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0082245-edf7-4e14-a5a2-a8150c05dec3",
   "metadata": {},
   "source": [
    "#### Define dS measure over the interface\n",
    "\n",
    "Interior facet integrals have no consistent ordering in `dolfinx`, resulting in the default assignment to `\"+\"` and `\"-\"` sides to be arbitrary. To restore a consistent ordering we define a custom `dS` measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7efcdb0a-85cd-4e7b-85f4-f4535793e623",
   "metadata": {},
   "outputs": [],
   "source": [
    "integration_entities_on_Gamma = dolfinx.fem.compute_integration_domains(\n",
    "    dolfinx.fem.IntegralType.interior_facet, mesh.topology, facets_Gamma)\n",
    "integration_entities_on_Gamma_reshaped = integration_entities_on_Gamma.reshape(-1, 4)\n",
    "connected_cells_to_Gamma = integration_entities_on_Gamma_reshaped[:, [0, 2]]\n",
    "subdomain_ordering = (\n",
    "    subdomains.values[connected_cells_to_Gamma[:, 0]] < subdomains.values[connected_cells_to_Gamma[:, 1]])\n",
    "if len(subdomain_ordering) > 0 and any(subdomain_ordering):\n",
    "    integration_entities_on_Gamma_reshaped[subdomain_ordering] = integration_entities_on_Gamma_reshaped[\n",
    "        subdomain_ordering][:, [2, 3, 0, 1]]\n",
    "integration_entities_on_Gamma = integration_entities_on_Gamma_reshaped.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb83bd0b-7d9c-4c2b-83f1-9ce85b5b762b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dS = ufl.Measure(\"dS\", domain=mesh, subdomain_data=[(2, np.array(integration_entities_on_Gamma, dtype=np.int32))])\n",
    "dS = dS(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff534460-ec4a-4393-9e53-4f52971fc48b",
   "metadata": {},
   "source": [
    "Check correctness of the subdomain measure by integrating a piecewise function defined as (1, 1) on subdomain 1, and (2, 2) and subdomain 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8818f6af-7587-4a8f-91ab-747e5e898987",
   "metadata": {},
   "outputs": [],
   "source": [
    "DG = dolfinx.fem.functionspace(mesh, (\"DG\", 0, (mesh.topology.dim, )))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55df1c39-179d-4704-87e5-ef448fd9ef35",
   "metadata": {},
   "outputs": [],
   "source": [
    "dg_function = dolfinx.fem.Function(DG)\n",
    "dg_function.interpolate(lambda x: np.full((mesh.topology.dim, x.shape[1]), 1.0, dtype=np.float64), cells0=cells_Omega1)\n",
    "dg_function.interpolate(lambda x: np.full((mesh.topology.dim, x.shape[1]), 2.0, dtype=np.float64), cells0=cells_Omega2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fd1ae96-34fa-4b78-849d-ab6a4b898275",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = ufl.FacetNormal(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6a808fb-3334-4566-b073-e2619e179762",
   "metadata": {},
   "outputs": [],
   "source": [
    "dg_function_check_from_left_subdomain = mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(dg_function(\"-\"), n(\"-\")) * dS)), op=mpi4py.MPI.SUM)\n",
    "dg_function_check_from_right_subdomain = mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(dg_function(\"+\"), n(\"+\")) * dS)), op=mpi4py.MPI.SUM)\n",
    "print(\"Check correctness when integrating from the left subdomain\", dg_function_check_from_left_subdomain)\n",
    "print(\"Check correctness when integrating from the right subdomain\", dg_function_check_from_right_subdomain)\n",
    "assert np.isclose(dg_function_check_from_left_subdomain, 6., atol=1.e-10)\n",
    "assert np.isclose(dg_function_check_from_right_subdomain, -12., atol=1.e-10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0292f97f",
   "metadata": {},
   "source": [
    "### With domain decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad602386",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function spaces\n",
    "V = dolfinx.fem.functionspace(mesh, (\"Lagrange\", 2))\n",
    "V1 = V.clone()\n",
    "V2 = V.clone()\n",
    "M = V.clone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2229e249",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define restrictions\n",
    "dofs_V1_Omega1 = dolfinx.fem.locate_dofs_topological(V1, subdomains.dim, cells_Omega1)\n",
    "dofs_V2_Omega2 = dolfinx.fem.locate_dofs_topological(V2, subdomains.dim, cells_Omega2)\n",
    "dofs_M_Gamma = dolfinx.fem.locate_dofs_topological(M, boundaries_and_interfaces.dim, facets_Gamma)\n",
    "restriction_V1_Omega1 = multiphenicsx.fem.DofMapRestriction(V1.dofmap, dofs_V1_Omega1)\n",
    "restriction_V2_Omega2 = multiphenicsx.fem.DofMapRestriction(V2.dofmap, dofs_V2_Omega2)\n",
    "restriction_M_Gamma = multiphenicsx.fem.DofMapRestriction(M.dofmap, dofs_M_Gamma)\n",
    "restriction = [restriction_V1_Omega1, restriction_V2_Omega2, restriction_M_Gamma]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76f66b48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define trial and test functions\n",
    "(u1, u2, l) = (ufl.TrialFunction(V1), ufl.TrialFunction(V2), ufl.TrialFunction(M))\n",
    "(v1, v2, m) = (ufl.TestFunction(V1), ufl.TestFunction(V2), ufl.TestFunction(M))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8ec8974",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define problem block forms\n",
    "zero = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0))\n",
    "a = [[ufl.inner(ufl.grad(u1), ufl.grad(v1)) * dx(1), None, ufl.inner(l(\"-\"), v1(\"-\")) * dS],\n",
    "     [None, ufl.inner(ufl.grad(u2), ufl.grad(v2)) * dx(2), - ufl.inner(l(\"+\"), v2(\"+\")) * dS],\n",
    "     [ufl.inner(u1(\"-\"), m(\"-\")) * dS, - ufl.inner(u2(\"+\"), m(\"+\")) * dS, None]]\n",
    "f = [ufl.inner(1, v1) * dx(1), ufl.inner(1, v2) * dx(2), ufl.inner(zero, m(\"-\")) * dS]\n",
    "a_cpp = dolfinx.fem.form(a)\n",
    "f_cpp = dolfinx.fem.form(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d88bf47",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define boundary conditions\n",
    "dofs_V1_partial_Omega = dolfinx.fem.locate_dofs_topological(\n",
    "    V1, boundaries_and_interfaces.dim, facets_partial_Omega)\n",
    "dofs_V2_partial_Omega = dolfinx.fem.locate_dofs_topological(\n",
    "    V2, boundaries_and_interfaces.dim, facets_partial_Omega)\n",
    "bc1 = dolfinx.fem.dirichletbc(zero, dofs_V1_partial_Omega, V1)\n",
    "bc2 = dolfinx.fem.dirichletbc(zero, dofs_V2_partial_Omega, V2)\n",
    "bcs = [bc1, bc2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a6958f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assemble the block linear system\n",
    "A = multiphenicsx.fem.petsc.assemble_matrix_block(a_cpp, bcs=bcs, restriction=(restriction, restriction))\n",
    "A.assemble()\n",
    "F = multiphenicsx.fem.petsc.assemble_vector_block(f_cpp, a_cpp, bcs=bcs, restriction=restriction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8990374e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve\n",
    "u1u2l = multiphenicsx.fem.petsc.create_vector_block(f_cpp, restriction=restriction)\n",
    "ksp = petsc4py.PETSc.KSP()\n",
    "ksp.create(mesh.comm)\n",
    "ksp.setOperators(A)\n",
    "ksp.setType(\"preonly\")\n",
    "ksp.getPC().setType(\"lu\")\n",
    "ksp.getPC().setFactorSolverType(\"mumps\")\n",
    "ksp.getPC().setFactorSetUpSolverType()\n",
    "ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=7, ival=4)\n",
    "ksp.setFromOptions()\n",
    "ksp.solve(F, u1u2l)\n",
    "u1u2l.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)\n",
    "ksp.destroy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88816060",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the block solution in components\n",
    "(u1, u2, l) = (dolfinx.fem.Function(V1), dolfinx.fem.Function(V2), dolfinx.fem.Function(M))\n",
    "with multiphenicsx.fem.petsc.BlockVecSubVectorWrapper(\n",
    "        u1u2l, [V1.dofmap, V2.dofmap, M.dofmap], restriction) as u1u2l_wrapper:\n",
    "    for u1u2l_wrapper_local, component in zip(u1u2l_wrapper, (u1, u2, l)):\n",
    "        with component.x.petsc_vec.localForm() as component_local:\n",
    "            component_local[:] = u1u2l_wrapper_local\n",
    "u1u2l.destroy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de3cd7ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(u1, \"u1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c7b65ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(u2, \"u2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9439d8c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(l, \"l\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e52653df",
   "metadata": {},
   "source": [
    "### Without domain decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e02cb9ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define trial and test functions\n",
    "u = ufl.TrialFunction(V)\n",
    "v = ufl.TestFunction(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef710299",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define problem forms\n",
    "a_ex = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx\n",
    "f_ex = ufl.inner(1, v) * dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5d85faf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define Dirichlet BC object on Gamma\n",
    "dofs_V_partial_Omega = dolfinx.fem.locate_dofs_topological(\n",
    "    V, boundaries_and_interfaces.dim, facets_partial_Omega)\n",
    "bc_ex = dolfinx.fem.dirichletbc(zero, dofs_V_partial_Omega, V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee02817d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve\n",
    "u_ex = dolfinx.fem.Function(V)\n",
    "problem_ex = dolfinx.fem.petsc.LinearProblem(\n",
    "    a_ex, f_ex, bcs=[bc_ex], u=u_ex,\n",
    "    petsc_options={\n",
    "        \"ksp_type\": \"preonly\", \"pc_type\": \"lu\", \"pc_factor_mat_solver_type\": \"mumps\",\n",
    "        \"mat_mumps_icntl_7\": 4\n",
    "    })\n",
    "problem_ex.solve()\n",
    "u_ex.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98d7eb35",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(u_ex, \"u\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "105ecd14",
   "metadata": {},
   "source": [
    "### Comparison and error computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fe214a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_ex1_norm = np.sqrt(mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex, u_ex) * dx(1))), op=mpi4py.MPI.SUM))\n",
    "u_ex2_norm = np.sqrt(mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex, u_ex) * dx(2))), op=mpi4py.MPI.SUM))\n",
    "err1_norm = np.sqrt(mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex - u1, u_ex - u1) * dx(1))), op=mpi4py.MPI.SUM))\n",
    "err2_norm = np.sqrt(mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_ex - u2, u_ex - u2) * dx(2))), op=mpi4py.MPI.SUM))\n",
    "print(\"Relative error on subdomain 1\", err1_norm / u_ex1_norm)\n",
    "print(\"Relative error on subdomain 2\", err2_norm / u_ex2_norm)\n",
    "assert np.isclose(err1_norm / u_ex1_norm, 0., atol=1.e-10)\n",
    "assert np.isclose(err2_norm / u_ex2_norm, 0., atol=1.e-10)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
